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Abstract. We study coupling functions that allow for persistent synchronization in 
,__, connected complex networks and any isolated system dynamics that possesses global 

C/3 solutions and bounded Jacobian evaluated along such solutions. We prove that the set 

^^ of coupling functions leading to stable synchronization is open and that any coupling 

(-h function whose linear part has eigenvalues with positive real part leads to system 

"£|2 towards synchronization. 
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1. Introduction 

Network synchronization is observed to occur in a broad range of applications in 
physics [TJ , neuroscience [21 El IH E] , ecology [6] , and life sciences [7j . During the last fifty 
years, empirical studies of real complex systems have led to a deep understanding of the 
structure of networks, the isolated dynamics of individual elements [8], the interaction 
structure [HI [10], and the interaction properties between oscillators, that is, the coupling 
function [HUE]. 

The stability of network synchronization is a balance between the isolated dynamics 
and the coupling function. Past research suggests that in networks of identical oscillators 
with interaction akin to diffusion, under mild conditions on the isolated dynamics, the 
coupling function dictates the synchronization properties of the network [131 HH [To] . 
However, it still remains an open problem to describe the class of coupling functions 
that lead the network to global synchronization for any isolated dynamics satisfying 
mild conditions. 

Our work contributes to the development a general theory for coupling functions 
that allow for persistent synchronization for a connected complex network. The coupling 
functions under consideration appear in a variety of synchronization models on networks 
(such as the Kuramoto models [12] and its extensions [El [20]). 

More precisely, we consider the dynamics of a network of n identical elements with 
interaction akin to diffusion, described by 

n 

x t = f(t,Xi) + a s ^2w ij H(x j -Xi) , (1) 

where a is the overall coupling strength, and the matrix W = (Wij)ij e {i ) .„ jn } describes 
the interaction structure of the network, i.e. Wy measures the strength of interaction 
between the nodes i and j. We make the following two assumptions for the function 
/:lxl ra 4 R m , and the coupling function H; R m -> R m . 

Assumption Al. The function f is continuous in the first argument and continuously 
differentiable in the second argument, and there exists an open set U C R m with C 1 
boundary that is e-inHowing invariant for some e > 0, i.e. the vector field f is pointing 
strictly inward at dll, uniformly bounded away from zero (see Definition^. Moreover, 
the Jacobian D x f with respect to x G R m is uniformly continuous and bounded on U, 
i.e. there exists B > such that 

\\DJ{t, x)\\<B for allteRandxeU. 

Note that if U is compact, then uniformity of the inflowing invariance condition as 
well as the uniform continuity of D x f and existence of a bound B follow automatically. 

Assumption A2. The coupling function H is continuously differentiable with 
H(0) = 0. We define T := DH(0) and denote the (complex) eigenvalues of T by 
Pi, i e {!,..., m}. 
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The network structure plays a central role for the synchronization properties. We 
consider the intensity of the i-th node Si = Y^j=i Wij, and define the positive definite 
matrix S := diag(Si, . . . , S n ). Then the so-called Laplacian reads as 

L = S-W. 

Let 6i, i G {1, . . . ,n}, denote the eigenvalues of L. Note that d\ = is an eigenvalue 
with eigenvector (1,1,..., 1). The multiplicity of this eigenvalue equals the number 
of connected components of the network. The last hypothesis is a crucial sufficient 
condition for synchronization. 

Assumption A3. We suppose that 

7 : = . min Re(0i&) >0, 

2<t<n, l<j<m 

where Ke(z) denotes the real part of a complex number z. 

The dynamics of such a diffusive model can be intricate. Indeed, even if the isolated 
dynamics possesses a globally stable fixed point, the diffusive coupling can lead to 
instability of the fixed point and the system can exhibit an oscillatory behavior |21j . 

Note that due to the diffusive nature of the coupling, if all oscillators start with the 
same initial condition, then the coupling term vanishes identically. This ensures that 
the globally synchronized state X\(t) — %z(f) — ■ ■ ■ — x n{t) — s {t) is an invariant state 
for all coupling strengths a and all choices of coupling functions H, and we call the set 

M :={xe (R m ) n :xi = ... = x n } 

the synchronization manifold. The main result of this paper is a proof that under the 
general conditions given above and a sufficiently large, the synchronization manifold S 
is exponentially stable. 

Theorem 1 (synchronization). Consider the network of diffusively coupled equations 
pi) satisfying A1-A3. Then there exists an a c = a c (f,T) such that for all coupling 
strengths 

a > — , 

7 

the network is locally uniformly synchronized, i.e. there exist C, 5 > such that if 
Xi(s) G U and ||xj(s) — Xj(s)|| < 5 for any i,j G {1, . . . ,n}, then 

\\xi(t) - Xj{t)\\ < Ce- (a7 - Qc)(t - s) ||xi(s) - Xj {s)\\ for all t > s . 

Hence, the characteristic relaxation time towards synchronization is qiven by — - — . 

Once the bound holds, the local dynamics will converge to the synchronization 
state and will be stable under small perturbations of the state. This means that the 
phenomenon of bubbling [22] and riddling [23] which leads to synchronization loss will 
not be observed, as opposed to the bounds that arise from the master stability function 
approach. Note that our approach is constructive and provides precise estimates for a c . 
These estimates depend on the properties of the matrices T and Df(s(t)). 
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Our second main result shows that this synchronization is persistent under 
perturbation of the isolated nodes. Thereto, consider a network of non-identical nodes 
described by 

n 

±i = fi{t, Xi) + a ^ W ij H i x j - Xi), (2) 

where fi(t,Xi) = f(xi) + gi(t,Xi). Note that in this case, the synchronization manifold 
S is no longer invariant. We will show in this paper that under the following additional 
conditions on the perturbation functions g^ i G {1, . . . , n}, the synchronization manifold 
is stable in the sense that orbits starting near the synchronization manifold M remain 
in a neighborhood of M. 

Theorem 2 (persistence). Consider a network of diffusively coupled equations and a 
choice of a > ct c jl as in Theorerrun Let M) describe a perturbation of the network such 
that 

11^(^,^)11 < £ g for allt el, x G W 1 and i6{l,...,n}. 

Then there exist e g ,5 > such that if \\xi(s) — Xj(s)\\ < 5 for any i,j G {1, . . . ,n}, 
then the network is locally approximately synchronized in the sense that the differences 
\\xi(t) — Xj(t)\\ still converge at an exponential rate towards close to zero. 



See page 16 for a more precise statement of this theorem. 



2. Discussion of the main results 

This section is devoted to relate our results to the state of the art, explain the 
assumptions and the ideas of the proofs. 

2.1. State of the art 

Recent efforts have focused on the role of the coupling function on the stability 
of the network synchronization. Pecora and collaborators have performed extensive 
experiments on various types of isolated dynamics and coupling functions [X6J to 
analyze the stability of synchronization. Pogromsky and Nijmeijer showed that if the 
linear coupling function is conjugated to a positive definite matrix, that is, symmetric 
and with positive eigenvalues, then it is always possible to synchronize a connected 
network [T7]. Ashwin and coworkers have introduce coupling functions that allow 
clusters of synchronization in the phase oscillators |18j . 

2.2. The assumptions 

Our assumptions are rather natural. Without assumption Al (existence of solutions) we 
could not speak of stability of the synchronized state. The second part of assumption Al 
(boundedness of the Jacobian) basically implies that with a finite value of a, choosing 
the coupling matrix properly we are able damp all the instabilities of the vector field and 
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obtain a stable synchronization state. Assumption A2 makes it possible to characterize 
the stability of synchronization behavior by the linearization of the H. Assumption A3 
guarantees that 7 is bounded away from zero. If this hypothesis is dropped, 7 may 
become negative and synchronization may no longer be possible. All these ingredients 
are used to prove the main result. 

2. 3. Ideas of the proofs 

To prove the result we map the synchronization problem to a corresponding fixed point 
problem. First, we start analyzing diagonalizable Laplacians. It is easy to show that if 
the coupling function is diagonalizable, then the diagonal dominance, see Proposition [5j 
implies that the synchronized state corresponds to a uniformly asymptotically stable 
fixed point. To obtain the claim for general coupling functions we make use of the 
roughness property associated with the equilibrium point, see Proposition |4j The main 
aspect here is to approximate the coupling function by a diagonalizable one while keeping 
control of the contraction rates. To conclude the proof for general Laplacians we prove 
that the set of diagonalizable Laplacians is dense in the space of Laplacians. From these 
results and the roughness property the main claim follows. 

2.4- The Assumption A3 

As far as we are aware condition A3 has not been put forward in the literature. 
Therefore, we would like to rephrase this condition when some further requirements 
are satisfied. 

The Spectrum of T is positive: If T has a spectrum consisting of only real 
eigenvalues then the condition A3 has a representation in terms of the Laplacian. In 
this case, the condition becomes 

Re(0i) > 

for every % ^ 1, as the Laplacian always has a zero eigenvalue. If the network is connected 
this eigenvalue is simple and in virtue of the disk theorem a sufficient condition for all 
eigenvalues to have positive real part is: 

The interaction strengths are positive, i.e., Wy > if i is connected to j, and zero 

otherwise. 

The Laplacian is Symmetric: This is the most studied case in the literature. 
The spectrum of the Laplacian is real. Assuming the network is connected we can order 
the eigenvalues as 

O = 0i<0 2 <03-"<0„ 

In this case the condition A 2 reduces to requiring the real part of the spectrum of T to 
be positive. Moreover, in this case the Laplacian L is diagonalizable by an orthogonal, 
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similarity transformation, which implies that the bounds on the persistence results 
reduces to 

\\xi(t) - Xj (t)\\ <Ck{Q) . 

«7 — a c 

as K 2 (P) = 1 and C is a constant depending only on the dimension of the isolated 
systems. 

3. Illustrations 

Before proving the main result, we consider two illustrative cases. First, two coupled 
nonautonomous linear equations. Then, we explore a case of coupled Lorenz systems. 

3.1. Nonautonomous Linear Equations 
Consider the nonautonomous linear equation 

Tt = A ^ x 



where 

-1-9 cos 2 (6t) + 12 sin(6t) cos(6t) 12 cos 2 (6t) + 9 sin(6t) cos(6t) 
-12 sin 2 (6t) + 9 sin(t) cos(6t) -1-9 sin 2 (6t) - 12 sin(6t) cos(6t) 

This is a classic example where the eigenvalues do not characterize the stability of 
the trivial solution of a nonautonomous equation. Indeed, the eigenvalues of A(t) are 
— 1 and —10, even independent of time, however, a direct computation shows that 

_ / e 2t (cos(6t) + 2sin(6t)) + 2 e - 13t (2cos(6£) - sin(6t)) 
X ^' ~ I e 2 '(cos(6t) - 2sin(6t)) + 2 e - 13 *(2cos(6t) - sin(6t)) 

is an unstable solution of the system. Hence, this system satisfies all the hypotheses. 
Consider now two diffusively coupled systems 

— - = A(t)xi + aT(x 2 - xi) (3) 

-^ = A(t)x 2 + «r(xi - x 2 ) (4) 

where T is a matrix. According to our main result, it is possible to synchronize these two 
systems for any coupling matrix with /3(T) > 0. Consider now the following coupling 
function 

(3 

Such a coupling function is in its Jordan form, and is non-diagonalizable. It is possible to 
synchronize these two systems by setting a properly. Consider the variable £ = x\ — x 2 , 
its evolution equation reads as 

| = [A(t) - 2oT]£. (5) 
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Our main result shows that the trivial solution of Eq. (p)]) is stable if a is large enough. 
We integrated Eq. (p)]) using a sixth order Runge-Kutta method with integration 
step 0.001. We computed the critical coupling a c as a function of (3, the result can 
be observed in Fig. [Tj The behavior of a c as a function of (3 appears to be intricate. 
For large j3 we obtain that a c tends to a constant, however, as we decrease (3 various 
changes of behavior can be observed. Although, the problem is linear, the critical 




Figure 1. We depict the critical coupling strength a as a function of A in a log — log 
scale. 



coupling strength depends nonlinearly on the parameter (3. 



3.2. Lorenz Oscillators 

Using the notation x = (u,v,w)*, the Lorentz vector field reads 

' a(v — u) ' 
f{x) = u(r — w)—v 
\ —bw + uv i 

where we choose the classical parameter values a = 10, r = 28, b = 8/3. The trajectories 
of the Lorenz eventually enter a compact set. Therefore, all trajectories of the system 
exist globally forward in time. Moreover, they accumulate in a neighborhood of a chaotic 
attractor |24j . 

Consider a network of three coupled Lorenz systems 

~dt 



(6) 



f(Xi) + a ^ W ij H ( X j ~ X i)i 
3=1 

the interaction matrix W is depicted in Fig. [2] 

We shall use two distinct nonlinear coupling functions, the first the associated 
matrix T is positive definite, whereas as in the second T is a Jordan block. The specific 
form of the coupling function can be seen in Fig. [3] 
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W 




Figure 2. The network and its weight matrix. The matrix L 
diagonalizable for every o^l, here, we choose a = 1/3. 



W is non- 



We integrated Eq. (|6j) using a sixth order Runge-Kutta method with integration 
step 10~ 4 , and computed the critical coupling a c as a function of (3, the result can 
be observed in Fig. [3] The behavior of a c is essentially different depending on the 
diagonalization properties of the associated I\ 



Jaap: I think that the result that a c ex /3 1 for /3 <C 1 is not surprising; essentially it shows that 
the coupling matrix with /3 = still induces synchronization. 
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Figure 3. Simulation results for the critical a c for the two coupling functions. For 
the first case, see left side, T = [31 is positive definite for /? > 0, and the behavior of 
a c does not depend significantly on f3. For the second case, T is a Jordan block with 
eigenvalues equal to /3. In this situation, for large values of /? the critical coupling 
a c appears independent of /3, as opposed to the small values of /?. In such case, the 
critical coupling scales as a c oc /3 _1 . 
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4. Notation and Background 

Here we quickly review results which are important to prove our Theorem. 

The object of our study is a network of identical systems with interactions akin to 
diffusion. We consider m-dimensional vector spaces IR m , and denote the elements x G M. m 
as column vectors x = (xi, 
space with the p-norms \\x\\ 



1 Xr 



(E 



where * stands for the transpose. We endow this 
:1 |xj| p ) . For simplicity, unless otherwise stated, 
the norm || ■ || represents the L^— norm ||x|| = maxj \xi\. Once we introduce a norm 
on the vector space R m , we also view the space of linear operators on M. m as a normed 
space equipped with the induced norm. The choice of norm in our case is immaterial 
as all our vector spaces have finite dimension. Furthermore, given a complex number z, 
we denote its real part as Re(z). 

4-1. Stability of trivial solutions 

Consider the linear differential equation 

dv rT , . 

* = mv (7) 

where v G M. m , for m > 1, and U(t) is a bounded and continuous matrix function. Recall 
that solutions of Eq. (J7l) can be written in terms of the evolution operator 

v (t) =T{t,s)v{s). 

The point x = is an equilibrium point. The time dependence in Eq. ([7]) introduces 
additional subtleties. Therefore, we want to have a stability condition that will also 
imply persistence under perturbations. Uniform asymptotic stability is such a condition; 
it is related to the evolution operator of the homogeneous equation in the following way. 

Definition 3. Let T(t, s) be the evolution operator associated with Eq. |?1). T(t, s) is 
said to be a uniform contraction if 

Vt>s:\\T{t,s)\\ <Ke-" (t - s) withr]>0. (8) 

The trivial solution of Eq. ([7]) is uniformly asymptotic stable if and only if the 
evolution operator is a uniform contraction. Moreover, uniform contractions have a 
rather important roughness property: they are not destroyed under small perturbations 
of the linear equations. 

Proposition 4. Suppose U(t) is a continuous matrix function on K and consider 
Eq. |?1). Assume that the fundamental matrix T(t, s) satisfies the exponential estimate 

Vt>s:\\T(t,s)\\<Ke p{t - s) (9) 

with K > and p G R. 

Consider a continuous matrix function V(t) satisfying 

su P \\V(t)\\=5 



Jaap: we 
should be 
careful since 
a change of 
norm adds 
to the overall 
constant; this 
can scale 
with n, m 
dimensions. 
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then the evolution operator T(t, s) of the perturbed equation 

f t = p{t) + V(t)}y, 
satisfies the exponential estimate 

Vt>s:\\f{t,s)\\ <Ke^- s \ 
where p = p + SK. 

Note that when p = —t] is negative, then T(t, s) is a uniform contraction, and 
setting 5 < r]/K we find as a corollary that p < 0, hence uniform contractions persist 
under small perturbations. 

This is a standard persistence result, and the proof can be found, for example 
in Ref. fI5\ Prop. 1]. There are various criteria to obtain conditions for uniformly 
asymptotic stability. We shall use the following criterion for diagonal dominant matrices 

Proposition 5. Let U(t) = \Uij(t)\ be a bounded, continuous matrix function on M. m 
and suppose there exists a constant rj > such that 

Re(U u (t)) + J2 \ U M\ <~V<0, (10) 

3=1, 

for all t 6 R and i = l,---,m. Then the evolution operator T(t,s) is a uniform 
contraction satisfying 

V£>s:||T(£,s)|| < e^^. 

The proof can be found in [25J. We use these fundamental results to prove our 
statements. 

5. Synchronization: Stability analysis 

Our goal is to prove that the synchronization manifold is locally attractive, and then 
to obtain bounds for the convergence of trajectories in an open neighborhood of the 
synchronization manifold towards the manifold. 

To this end, we first obtain equations that govern the dynamics near the 
synchronization manifold. The dependence of / on time is not crucial in the proof, 
hence omitted for simplicity. Recall that we defined F := DH(0). Using a tensor 
representation we can write the n equations as a single equation. First define 

X = col(xi,...,x n ), 

where col denotes the vectorization formed by stacking the column vectors Xj into a 
single column vector. Similarly, we define 

F(X)=col(f(x 1 ),---J(x n )). 

We can analyze small perturbations away from the synchronization manifold in terms 
of the tensor representation 

x = i® s + e, (ii) 
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where ® is the tensor product, and 1 6l" such that 1 = col(l, • - - , 1), which is the 
eigenvector of L corresponding to the eigenvalue zero. Note that 1 <g) s defines the 
synchronization manifold, and we view £ as a perturbation to the synchronized state. 



Let us briefly expand a bit on the coordinate splitting in (11). Our state space 
(R n (g) R m , ||-||) can be canonically identified with IR nm , which we will use for shorter 
notation. The norm need not come from an (Euclidean) inner product. For example, it 
can be the maximum over the Euclidean norm \\-\\ E on each of the R m spaces, i.e. 

\\(xi, . . . ,x n )\\ = max \\xi\\ E where x, G M m . (12) 

l<j<71 



The coordinate splitting (11) is associated to a splitting of lR nm as the direct sum 
of subspaces 

R nm = M®N (13) 

with associated projections 

tt m : R nm -> M, tt n : R nm ->• N. 

The subspaces M, iV C R nm are determined by embeddings from IR m and IR*™ -1 )" 1 
respectively, induced by the Laplacian matrix L on 1". 

Let us for the moment use the simplifying assumption that L is diagonalizable with 
eigenvectors I,v 2 , ...,v n . Then the subspaces M,N have natural representations in 
terms of these eigenvectors as 

M = span(l) ® R m , N = span(w 2 , ...,v n )<8> K m . 

This means that we have 'natural' embeddings that induce coordinates on these 
subspaces: 

i M :K m -^M, si-H®s = col(s,...,s), 

n 

If we drop the assumption that L is diagonalizable then we lose the natural choice of 
an embedding for N. 

The norm ||-|| on R nm can be restricted to the subspaces M, N and induces norms 
on the 'coordinate' spaces R m , R( n-1 ) m by pullback under the embeddings. That is, the 
induced norm on IR m is given for example by 

IMIr™ = lkw(*)ll = II 1 ® s ll> 

which is equal to \\s\\ E if ||-|| is chosen to be the maximum of the Euclidean norms as 



in (12). Henceforth we shall identify s E R m with 1 <g> s E M under this isometry lm- 

In the next proposition we represent the dynamics in terms of s E M and £ E N. 
The idea is that X is in an open neighborhood of M, i.e. that £ is small. 



Jaap: K 
conflicts 
with use as 
constant in 
exponential 
estimates 
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Proposition 6. The perturbation £ satisfies the equation 

f t =K(s)Z + R(s,0, (14) 

where 

K(s) = I n ®Df(s)-a(L®r) 

and R(s, £) is the remainder satisfying the following property: for any e > 0, there 
is a 5 > such that for all ||£|| < 5 one has ||i?(s,£)|| < £||£|| uniformly in s, 
i.e. R(s,£) G o(||£||) uniformly in s. 

Note that this equation for £ still depends on s; later we shall provide stability 
bounds for it, independent of s. 

Proof. Since H is smooth, we have by Taylor's theorem 

H(x) = T x + r(x) 
with ||r(a;)|| < ellxll for \\x\\ < 5. Now we define 



R H {X)i = Y^ w ij r ( x i ~ x j) 



i=i 



J2 Wijrfai 1 ®s + 0-Pj(1L®s + 0) 



n 

= J2 w ^<p^-pA0), 

where Pi denotes projection onto the ith R m tuple in R nm . So we see that Rh(X) = 
Rh{0 does not depend on s G M, and it satisfies the estimate 

n 

\\Rh(Z)\\ < max (X;|Wi i |)e2||e|| when ||£|| < 8/2. 
i=i 
Recalling that L^ = SijSi — W^, the coupling term can then be rewritten as 

n n 

Y, W l3 H{x, - Xi ) = -^2 L tJ T Xj + R H (0i = (L <g> F)X + R H (£). 

3=1 3=1 

Let us now look at the term F(X) describing the dynamics of the uncoupled nodes. 
We Taylor expand X = 1 ® s + £ around 1 (g) s and find 

F(l <g> s + = F(l <g> s) + £>F(1 <g> s)£ + i? F (s, 

= l®/(s) + J n ® £>/(«)£ + ^(5,0, 

and we use a mean value theorem estimate to obtain that ||i?p(s, £)|| < e||£|| when 
||£|| < 5 and Df is uniformly continuous. 
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Thus we recover the full differential equation for the system, rewritten in the 
coordinates (s, £) G M © iV as 

d aW = 1 ®Jt + < ft =1 ® /(S) + In ® jD/(s)e ~ a{L ® r)e 

+ R F (s,0 + aR H (0- (15) 

Next, we project both sides onto the spaces M, N to recover the separate differential 
equations for s, £ respectively. This leads to 

s = f(s) + n M (R F (s, + aR H (0), (16) 

i = K(s)£ + MRf(s, + «i? ff (0), (17) 

where 

tf(s) = J n <g>D/(s)-a(L<g>r). 

Note that both J n ® Df(s) and L <g> T preserve the subspaces M and iV, since J n and L 
preserve both span(l) and span(v 2 , . . . ,v n ), thus the projections can be dropped there. 
Furthermore L t = means that the term (L <g> T)(l ® s) disappears from the equation 
as well. □ 

Before we can continue to analyze the flow for £, we need to control s. 
Assumption yO that the single node system has a uniformly inflowing invariant set 
U C M m leads to a similar result close to the synchronization manifold M in the coupled 
network. This guarantees that the solution for s stays inside U C M m , irrespectively of 
the precise details of its flow and the dependence of the flow on £. 

To prove this result let us first introduce some additional notation. The following 
definition of uniform inflowing invariance was already referenced in assumption Al. 

Definition 7 (^-inflowing invariance). Let x — f(x) describe a dynamical system with 
x G lR m and let U C 1R" 1 have C 1 boundary. Then we call U an e-inflowing invariant 
set for f if we have at each point x G dll with inward-pointing normal vector n x that 

(n x ,f(x))>e. (18) 

Secondly, we need the concept of a tubular neighborhood of a submanifold. 

Definition 8 (^-tubular neighborhood). Let t <g) U C W 1 " 1 be a subset of the 
synchronization manifold with C 1 boundary and let r\ > 0. Then we call 

Ur, = {i®s + S\seU,SeN, U\\<v} (19) 

an rj-tubular neighborhood of t <g) U. 

The following lemma implies that if a solution curve (s(t),£(t)) leaves U v , then it 
must do so by ||£(£)|| growing larger than r]. As a result we can ignore the dependence 
on s of the flow for £ as long as we have ||£|| < r\: the coupling does not depend on s, 
and by assumption we have for s G U that the terms Df(s) and Rp(s,£) are uniformly 
bounded and small, respectively. This lemma is formulated with a general perturbation 
G to allow its reuse in the persistence proof. 
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Lemma 9. Let the system x = f(t, x) satisfy assumption Al with associated e-inflowing 
invariant set U C IR m . Let X = F(t,X) describe the dynamics of n uncoupled copies of 
this system and let G:lx W im — y IR nm be a perturbation to F such that for some r > 
and 5 > it holds that 

sup \\G(t,X)\\<5<e/\\7c M \\. 

Then there exists anO < rj < r such that solution curves (s(£), £(£)) of the system defined 
by F + G can only leave the tubular neighborhood U v through 

d C y l u v = {i®s + z\seU, U\\=v}- 

Proof. Let us consider some < rj < r to be fixed later. Since U v is cylinder-like, its 
boundary consists of two parts: 

dU v = dcyiUrf U dsideU^, 

that is, the points 1 £g> s + £ where ||£|| = n and those where s G dU. We only need to 
consider the second case, that is ^ide^- 

If n is an inward pointing normal vector at s G dU, then 1 <g> n points inwards to 
dsideUn at 1 <g> s + £. Note that 1 <g> n cannot be said to be normal to d S id e U v since a 
natural inner product on IR nm is missing. Instead we show that the component of F + G 
projected along M has positive inner product with n, and thus F + G points inwards 
at dU v . Here we use the isometry lm to endow M with the inner product of M. m . 

From the bounds \\DF\\ = \\Df\\ < B and ||G|| < S it follows that 

(n, tt m [F(1 (8) s + + G(l (8) 5 + 0]} 
= (n, f(s)) + (n, n M [DF{l ® s + r £)£ + G(l ® s + 0]> 
>e- 11^11(11^1111^11 + ||G||) 
>£-|Km||(S77 + «J). 

We applied the mean value theorem with r G (0, 1) as interpolation variable. Since 
IKm||^ < £ it follows that there exists an rj > such that F + G points inwards 
everywhere at d S id e U v . □ 



To apply Lemma M to our system, we note that the terms in ( 15 ) that do not arise 
from F are all linear or higher order in £, so if we denote these collectively by G. Then 
G can be made smaller than a given S > by restricting to a sufficiently small tubular 
neighborhood U r for some r > 0. Thus, if we ensure that ||£(£)|| < r\ for all t, for a 
solution curve that starts inside U v , then it stays inside U v . 

5.1. Diagonalizable Case 

For clarity we first tackle the case where L and T are diagonalizable. We obtain the 
general case via perturbations and the roughness of dichotomies. 
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Lemma 10 (Diagonalizable case). Let the network satisfy the assumptions A1-A3. 
Moreover, assume that the network is connected and that the Laplacian L and linearized 
coupling T are diagonalizable: 

L = PAP' 1 and T = QGQ' 1 

where A = diag(8i,8 2 , . . . , 9 n ) and G = diag(0i, . . . , (3 m ) are the eigenvalue matrices of 
A and F respectively. Let s denote a solution of s — f(t, s) with s(£o) = s o G U and let 
K = K(s(t)) denote the associated linear flow for £ E N as defined in Prop. [6l Then 
there exists an a c such that for all 

a c 
a > — 

7 
the linear evolution operator T(t, t ) = T(£, £o)|at for £ is a uniform contraction with 
estimate 

Vt>t : \\T{t, t ) || < Cn(P)e~ {a ^- ac)it - tQ) (20) 

where C = C(Q) > 1. 

Proof. Note that S = P <S> Q is an invertible matrix that diagonalizes L £g> T and the 
change of coordinates 

K = S- 1 KS = I n ® Q- 1 Df(t, s(t)) Q-aA®G (21) 

reduces K to m-block diagonal form. Thus we have 

iT = 0^ where K t = Q' 1 Df(t,s(t)) Q - aOiG. 

8=1 

Since K is block diagonal, it preserves the splitting IR nm = ®" =1 IR m , and hence its 
associated flow T(t, t ) will also be of the form 

n 

f(t,t ) = ($f l (t,t ), (22) 

where each Tj(t,to) is the flow of Ki. Note that the spaces M m correspond to Vi ® M. m 
under the coordinate change S and in particular M = S M.™ where M™ denotes the first 
copy of the n spaces lR m . Thus, restricting K to N corresponds to restricting K to the 
blocks i > 2. 

Now for each block we obtain 

—7- = ki(s)yi = [J -a 9i G]yi 

where J := Q^ 1 Df(t, s(t)) Q. By assumption Al we have that s(t) G U for all t > to, 
and also 

sup \\Q- l Df{t,s)Q\\ <k{Q)B. 
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Now, to apply Proposition [5] to this equation, we want 

m 

Re[J fcfc - aBiGkk] + Yl l J *i(*)l < ° 



3=1, 



to hold. Since Re(Jkk) < \Jkk\ we therefore find that 



a > 



Re(6' i G fcfc ) 

Noting that Y^T=i l"A?l — K (Q)B and 7 < Re(#jGfcfc) we obtain as sufficient condition 
that a > a c := (k(Q)B)/j. Therefore, the evolution operator Xj(£,£ ) satisfies 

\\fi(t,t )\\ < (7 e -(«T-«(Q)B)(t-to) 

where C = C(m) is a constant depending only on the dimension of the isolated nodes. 



Finally, using Eq. (22) and changing back to the original coordinates it follows that 

|s( , 0f i (Mo))s , - 1 l. 



lim«)kl 



< ft(S') max Tj(t, ti 



j>2 



o, 



< k{P) k{Q) C{m) e -(^-K(Q)B)(t-t ) _ ( 23 ) 

Note that 5* _1 maps M and A" onto the first and last n — 1 of the m-tuples in IR nm 
respectively, hence the restriction to iV reduces to a direct sum over i > 2 after 
conjugation with 5*. Also, we are only considering the parts of S and S~ l restricted to 
S~ 1 N and A" respectively, but we can simply estimate k(S\s-i n ) < k(S). □ 



Jaap: prove statements of Theorem [l] now: that U is uniformly attracting (including nonlinearities) 
and that the exponential rate is a c — aj < 0. 



6. Proof: A result on Persistence 

Recall that we consider a network of almost-identical nodes described by Eq. (J2l): 



Xi == h(t,Xi) + a > Wan (Xj - Xi), 



3=1 

where fi(t,Xi) = f(t,Xi) + gi(t,Xi). We shall now prove the following more precisely 
restated version of Theorem [2] 

Theorem ^f (persistence). Consider a network of diffusively coupled equations and a 
choice of a > ct c ll os in TheoremUl Let M) describe a perturbation of the network such 
that 

\\gi(t,x)\\<£ g for allt 6R, x G M m and i G {!,..., n}. (24) 
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Then there exist C, 5 > and Eq > such that if e g < Eq and \\Xi{s) 
any i,j G {1, 



xAs 



. . ,n}, then 

\ Xi (t) - Xj {t)\\ < Ce- {a ~<~ ac){t - s) \\xi{s) - Xj {s) 



+ 



Ce r , 



< 5 for 

(25) 



a c — oey 
for all t > s. 

Note that the proof does not specifically depend on the fact that the perturbations 
gi to the nodes are decoupled; the function G below can depend in an arbitrary way on 
the total state X as long as it is uniformly small. 

Proof. Denote by 

G(t,X) = col{gi{t,xi),...,gi(t,x n )) 

the perturbation for the whole system. Note that ||G|| < e g holds. When IR n is endowed 
with the maximum norm, we have the following equivalence of norms: 



FAT 



I I^H sup ||x 4 -xj< HfH 



l<i,j<n 

\X±, ■ ■ ■ , x r 



(26) 



for any X = l®s + ^ = (xi, . . . ,x n ) G M. nm . Thus, the estimates ||xj(s) — Xj-(s) || < 5o 
in Theorem [2] imply ||£|| < ||7rjv||^o- Since the assumptions A1-A3 are uniform in time, 
we can again suppress explicit time dependence and restrict to the case s = 0. 



The addition of G modifies the differential equations (16) and (17) for s and 



£. Lemma M guarantees that there exists again an 77-tubular neighborhood U v such 
that solutions of the complete system for (s,£) cannot 'escape along s\ when e g ,rj are 
sufficiently small. 

The differential equation for £ is given by 

£ = K{s)£ + R(s, + ir N (G(l <g> s + £)), (27) 

c.f. Proposition [6J We denote by e(S) the continuity modulus of R. 



Variation of constants yields that solutions of (27) satisfy 



at) = T(t,0)m+ / T(t,r) i2( S (r),£(r)) + G(l®8 + £) 



dr 



with s(t) a solution of the modified equation for s and T(t,r) the evolution operator 
associated to K(s(t)). 

Let us now assume that ||£(0)|| < 1 1 7Ttv 1 1 <^o and ||£(£)|| < 8\ < r\ for all t > 0. From 
the proof of Theorem hi it follows that \\T(t, r)\\ < Ce f> ^ T " > with p = a c — ccy < 0. Then 
we readily estimate 

-t 



\\m\\<Ce» t \\'K N \\5 Q + / Ce^- T \e(S 1 )U( 



\-K N \\s q ) dr 



(28) 




C 



< Ce pt \\ir N \\6 + —(e(S 1 )S 1 \\7i N \\e g ). 



□ 



Jaap: do an explicit Gronwall-type estimate for (28) to recover (251 
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7. Normal hyperbolicity 

In the previous section we proved persistence of synchronization under the condition 

a c — ccy < 0. 

The left-hand side of this inequality is precisely the rate of exponential attraction 
towards the synchronization manifold. Indeed, this condition persists under small 
perturbations, so solutions of the perturbed system still converge to approximately the 
synchronization manifold. However, there need not exist an invariant manifold in the 
perturbed system that precisely attracts all nearly solutions. For this to be true, the 
original synchronization manifold must be a normally hyperbolic invariant manifold, 
see J27J EH] for fundamental definitions and results and e.g. [2H] for normal hyperbolicity 
in synchronization of networks. 

The additional condition for normal hyperbolicity to hold is that the flow 
contracts less along the synchronization manifold than that it contracts in the normal 
directions. The estimate a c on the Jacobian of / was used previously to estimate 
the maximum expansion that the flow of / can detract from the contraction in the 
normal direction; similarly — a c provides an estimate for the maximum contraction along 
the synchronization manifold. Thus the so-called spectral gap condition for r-normal 
hyperbolicity is satisfied when 

a c — aj < —r a c with r > 1. 

When the synchronization manifold satisfies this normal hyperbolicity condition, then 
the stronger result holds that it persists under arbitrary C l small perturbations. This 
means that not only solution curves converge to close to the original synchronization 
manifold, but that there exists a persistent C r synchronization manifold 

M = {xi = hi(s), s6K ra ,l<i<n} 

where the functions hi are uniformly C 1 close to the identity function, i.e. M is close 
to M. Solution curves of the perturbed system converge at an exponential rate p close 
to a c — aj to this manifold M, and moreover a stronger 'shadowing' or 'isochrony' 
property holds that any solution curve X(t) that converges to M, actually converges at 
exponential rate p to a unique solution curve X^(t) on M in the sense that there exists 
a C such that for all t > 

\\X{t)-X a {t)\\<Cef*. 

Sharper estimates can be obtained however, when it is known that the 
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